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DYNASTY  U,  A  NONLINEAR  IMPUCIT  CODE  FOR 
RELATTVISnC  ELECTRON  BEAM  TRACKING  STUDIES 


I.  Introduction 

Computational  models  that  are  linearized  about  an  axlsymmetrlc  equilibrium 

have  proven  Invaluable  In  the  study  of  Instabilities  of  propagating  beams.  For 

example,  linearized  (monopole/dlpole)  particle  simulation  codes^  as  well  as 

2 

linearized  codes  with  simplified  models  of  beam  dynamics  have  been  used  for  hose 
Instability  analysis.  However  there  are  other  areas  of  beam  propagation 
phenomenology  In  which  large  departure  from  axlsymmetry  are  essential  and 
Inevitable,  and  In  which  such  linearized  models  can  be  quite  misleading.  For 
example.  In  the  presence  of  transverse  external  magnetic  or  electric  fields, 
different  "slices"  of  the  beam  (defined  by  the  distance  r,  from  the  beam  head)  may 
be  subject  to  different  accelerations  F/Corf),  where  F  Is  the  force  and  Yt  vhich 
may  be  a  function  of  Is  the  relativistic  factor.  In  such  situations,  the  beam 
may  tear  apart,  or  alternatively  it  may  reach  a  non-axlsymmetric  equilibrium  in 
which  Its  cohesive  self-forces  counterbalance  the  sheared  external  forces.  As 
another  example,  we  have  recently  shown  that  there  are  Important  situations  In 
which  the  beam  head  Is  In  a  grossly  unstable  equilibrium  (unstable  to  simple 
zero-frequency  transverse  displacements)  when  It  Is  propagating  on-axls  In  a 
density  or  conductivity  channel,  but  In  which  another  and  more  stable  equilibrium 
exists  with  the  beam  head  off-axis  by  a  prescribed  amount.  Thus  axlsymmetry  Is 
broken  and  the  beam  tends  to  track  the  channel  by  riding  down  the  channel  walls, 
rather  than  the  channel  center.  Moreover,  the  equilibrium  displacement  from  the 
channel  axis  Increases  from  the  head  of  the  beam  to  Its  tall,  Increasing  the 
departure  from  axlsymmetry. 

To  treat  problems  such  as  these.  It  Is  necessary  to  develop  a  fully  three- 

dimensional  solution  of  the  electromagnetic  field  equations;  monopole/dlpole 
Mamuscript  approved  May  25,  1983. 


approximations  will  not  do.  In  dealing  with  very  highly  relativistic  paraxial 

beams  we  shall  assume,  however,  that  it  is  sufficient  to  solve  the  three- 

dimensional  form  of  Lee  s  reduced  field  equations  ,  rather  then  the  complete 

Maxwell  equations.  This  paper  deals  mainly  with  the  numerical  methods  which  we 

have  developed  to  solve  the  field  equations.  To  date,  we  have  only  addressed  the 

problem  of  determining  non-axlsymmetrlc  beam  equilibria,  in  which  the  beam  Is 
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specified  by  a  prescribed  density  profile  n^(r-Y(i;),  5).  Each  slice  of  the  beam 
Is  allowed  transverse  (non-lnf Inltesimal)  displacement  Y(;),  but  Is  not  allowed 
to  distort.  No  particle  dynamics  and  no  time  dependences  have  been  studied,  to 
date.  Eventually  we  shall  have  to  study  the  stability  of  these  equilibria  to 
Internal  distortions  and  to  dynamic  modes  such  as  hose,  but  this  can  only  be  done 
after  the  correct  equilibria  have  been  determined,  a  task  which  Is  by  no  means 
trivial  either  computationally,  mathematically  or  physically. 

We  have  already  discussed  the  code  Dynasty-I,  which  was  developed  as  our 
first  cut  to  study  the  Interaction  between  the  beam  and  the  channel.  It  solves 
the  full  set  of  Lee  field  equations  without  linearization.  The  radial  dependence 
Is  solved  by  finite  difference  with  nonuniform  grids.  The  theta  dependence  Is 
solved  by  fast  Fourier  transform.  The  advantage  of  this  method  Is  that  the  theta 
derivatives  are  more  than  second  order  accurate^.  Unfortunately,  Dynasty-I  Is 
severely  time  step  limited  because  It  uses  an  explicit  scheme.  This  can  be  seen 
from  the  following:  the  Lee  field  equations  are 


where  n  =  o  •  Both  n  and  a  behave  like  diffusion  coefficients  In  the 
equations.  Near  the  head  of  the  beam  n  Is  very  large,  but  near  the  tall  a  Is 
very  large.  The  difference  can  easily  be  three  orders  of  magnitude.  So  any 
explicit  scheme  that  takes  care  of  the  solution  at  one  end  will  have  severe  time 
step  limitations  at  the  other  end. 

Recently  we  have  developed  a  new  code  called  Dynasty  II.  It  Is  similar  to 
Dynasty  I  except  It  solves  the  field  equations  Implicitly  and  so  the  time  step 
limitation  Is  removed.  This  code  Involves  the  solution  of  a  bl-trldlagonal 
system  with  complex  matrix  coefficients.  The  complexity  In  calculation  Is  offset 
by  the  fact  that  for  the  first  time,  one  Is  able  to  find  the  self  consistent 
fields  and  the  force  acting  on  the  beam  from  head  to  tall  for  arbitrary  beam 
displacement . 

We  will  not  present  the  physics  results  Investigated  by  Dynasty  II  here; 
they  will  be  reported  In  later  publications.  Instead,  we  shall  concentrate  on 
the  numerical  formulation  of  the  code.  In  section  II,  we  review  briefly  the 
convolution  sum,  and  In  Sec.  Ill  we  apply  It  to  a  simple  diffusion  equation.  In 
Sec.  IV,  the  method  of  solving  the  field  equations  Is  described.  In  Sec.  V,  we 
discuss  some  applications. 


II.  Review  of  Convolution  Coding 

Let  A(9)  and  B(6)  be  two  functions  that  depend  on  9.  The  expansion  of  the 
product  of  A  and  B  In  a  finite  Fourier  series  Is 


A(9)  B(9)  -  5:  A  e"^^®  I  B 

1— N  *  m— N 


N  N 


I  I  A  B 


-1(£  +  m)9 


£— N  m— N 


Let  n  »  £  +  m,  then  (1)  becomes 


N  m+N 

y  y  A  B  e' 

m“-N  n“m-N  ^ 


We  want  to  convert  (2)  to  a  form  like  ^  (  )  e  .  Breaking  the  double  sum  to 

o*-N 

two  sums  and  Interchanging  the  direction  of  Integration,  we  obtain 


2N  N 


A  B  -  (f  y  y  ]  +  [  y  y  ii  a  b 

n“-2N  m*-N  n"0  m*n-N  ™ 


Equation  (3)  can  be  approximated  by  dropping  the  higher  order  modes,  l.e.  we 
2N  N 

replace  y  by  y  .  This  can  be  justified,  provided  that  the  magnitude  of  a 
n— 2N  n— N 

transform  goes  down  rapidly  with  higher  mode  number.  Equation  (3)  becomes 


N  N 


A(9)  8(9)  -  {[  y  I  ]+[  I  I  ]lAn_^B^e 
n—N  m*-N  n*0  m“n-N 
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III.  Implicit  Formulation  of  a  Diffusion  Equation  by  the  Transfoinn  Method 
The  RHS  of  a  diffusion  equation. 


H  (r,9)  -  n(r,8)  LtA(r,9)] 


(5) 


with  L  =  ,  can  be  expressed  as 


0  n+N  N  N 

“S  ■  {[  I  I  1  *  [  I  I  ]1 

n— N  m"-N  n"0  m*n“N 


-in9 


(6) 


In  matrix  form,  (6)  becomes  (dropping  e 


+  n_N(LA)^ 


+  0 


ni(LA)_jj  +  %^^^^-N+l 

#  • 

ni,(LA)_jj  +  nj,.i(LA).j^^+ 


+  •  •  • 


^(LA)o  +  •  •  •  +  nQ(LA)jj 


(7) 


A  typical  (I-A)|jj  term  can  be  expanded  to 


(8) 


Using  a  nonuniform  grid  in  radial  direction,  i.e.  using  uniform  spacing  in  y 


ri_i 


n.!,) 


-N 

Sn+1 


o  J 


1+1 


•  *  ^-N^i,o^ 


-N 


-M+1 


[;  J 


n  „) 


-N 

Sn+i| 


o 


1-1 


(15) 


Similarly  we  can  rewrite  all  the  rows  In  (7)  In  a  form  such  as  (15). 
Consequently,  the  RHS  of  (5)  can  be  expressed 


m 


!•! 


^1+1-^ 


A  +  T 
Af  -r 


1-1 


(16) 


where  A^_j^  are  column  matrices  with  the  first  element  equal  to  A 

the  last  element  equal  to  Am,  and  R, ,  S . ,  T.  are  banded  matrices  with 


•  • 


n, 


-l  *  •  •  ^_N 


^2  %  *  •  * 


n^.i  •  •  •  n.^ 


\*  *  •  \  \ 


and  aj^  and  are  diagonal  matrices  with  their  elements  equal  to  and 
y^  respectively.  0^^  is  also  a  diagonal  matrix,  but  each  element  is  different 


Going  back  to  (5),  we  finally  obtain 


^i  h  1  [h  ^1+1  ^  ^i  ^i  ^  ^i  Ai_iJ 
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2  Li  i+1 


i  i 


i-lj 


5  = 


n+1  n 
C  -  C 


(22) 


The  scheme  described  in  (21)  is  similar  to  a  Crank-Nicholson  method  and 
2  2 

is  0(6)  and  0(A)  •  Now  we  are  ready  to  tackle  the  Lee  field  equations,  as 
discussed  in  the  next  section. 


Implicit  Formulation  of  the  Lee  Field  Equations  by  Fourier  Transform 
We  can  cast  Lee's  field  equations^  in  to  the  following  form: 


where  n 
system. 


where 


2 

—  -  n  (A  +  <^)  +  n  (23) 

2  SA  O 

(24) 

=  a  Our  aim  here  is  to  reduce  (23)  and  (24)  to  a  bi-tridiagonal 
The  conductivity  equation  in  normalized  unit  is 

-  V  T  4.  V  2 

H  o  -^2?  o  , 


-  1.4653 

-  1.7  X  10“^ 

APS^ 

'^i  *  2  ’ 

1  +  BS  +  CS  +  DS 
S  »  E^/P^, 

A  -  1.423  X  10'^, 

B  -  9.179  X  10~^, 


■ 


2.656  X  10 


-10 


D  =•  2.820  X  10 


E  and  P  are  the  local  electric  field  ?nd  pressure  respectively.  The  equation  for 

0  or  n  is  not  solved  implicitly  because  it  includes  nonlinear  terms  due  to 

avalanche  and  recombination,  but  we  can  always  time-center  the  o  equation  to 

2 

ensure  accuracy  to  0(6)  .  For  simplicity,  let  us  define 


a 


- 

""  3r 


«• 


a 


-  l£ 

■  39 


(25) 

(26) 


(24)  can  be  rewritten  as 


We  shall  now  transform  (23)  and  (27)  teirm  by  term.  From  (23),  we  have 


(30) 


*1  t®l 


n+1  /2 

The  term  (n  J,  }  is  calculated  explicitly.  From  (27),  we  have 


2  8A  .-1  r  r-  rn+1  -  rn+1  -  rn+1 

^  i  tVi+1  hh  ^iVi. 


-[\  \+i  + 


1^"  (8yiA)“^  ~a[  (ii+r  ii-i)''].  (32) 


a"  3  ^  i  -  r-n+1  ^  -Til 

T  Te  *  TT  ‘^i  “  t^i  •"  ^il‘ 

r  2y^ 

-  1/2 

where  i  =  (-1)  to  avoid  confusion  with  index  i.  Also 


(j)  2  L^i  *1*1+1  ^  ^i  ^i  ^  1  ^i-lJ 


In  the  above,  m  is  a  diagonal  matrix  with  element  runs  from  -N  to  N,  and  R^, 
and  are  defined  in  (17)-(19)  with  replaced  by  a^.  Putting  (28)-(30)  in 
(23)  and  (31)-(34)  in  (27),  we  get  the  desired  bi-tridiagonal  form: 
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(48) 


(4)  ^  ^1 

-  -  -57:7  -T-* 


4.fT  +  _ i.)  a”  +  _ i.  A  " 

2  ^i+1  ^  ^  2-  ■*  ^  2  ^i- 


1 
6  S 


.  6  T, 

1-1  .n  .  i  .  n 


+  1  (Rt^i+l+  \  Cl)  +  (nJt)""''^6 


>  n+l / 2 


(49) 


5  a.  6R 


/0\  _  \J  W* 

■*1  ■  »1  ^+1*  h  ^-r  (TTi * -r) 


i^  -rn 
i+1 


7x  "  4  S  6o.  5T 

-  I  4  ®i  “  "  2  )  ‘•’i  "  Uy.A  "  2  ) 

2y.  •'i 


(50) 


where  I  is  an  unit  diagonal  matrix.  The  Inversion  method  of  the  bl-trldlagonal 
system  Is  given  In  the  appendix. 

We  use  conducting  boundary  conditions  at  the  outer  boundary  l.e. 


A„(R)  -  0, 

V(R)  -  0, 
m 


(51) 


for  all  m.  Near  the  origin,  we  require  both  the  electric  and  magnetic  fields  be 
finite  and  continuous.  Noting  that  for  small  r. 


A™(r).  ^„(r)  «  r“,  (52) 

ni  ts 

we  have  at  once 

A„(0)  -  0  -  1)„(0),  for  m  >  1,  (53) 

m  m 

3A  (0)  3A  (0) 

-^-0--^,  form-0.  (54) 
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V.  Applications 

To  date.  Dynasty  II  has  been  used  only  to  calculate  beam  equilibria,  using  a 
simple  envelope  model  of  the  beam.  The  beam  current  density  Is  specified  to  be 
of  the  form 

fir  -  V(5)ll 

^  [  ~V(T)  J- 

where  f  Is  a  specified  radial  profile  shape  (usually  Bennett),  a(;)  Is  the  beam 
radius,  and  Y(;)  Is  the  transverse  displacement  of  the  slice.  In  most  cases,  we 
have  specified  a(;)  and  used  the  code.  Iterating  back  In  ;  from  the  beam  head,  to 
calculate  the  equilibrium  value  of  Y(;).  We  have  already  mentioned  briefly  some 
applications  of  the  code  In  Sec.  I.  We  shall  elaborate  on  these  areas  In  some 
detail  here. 

Previous  work^  has  shown  that  In  the  electrostatic  regime  at  the  beam  head 

If  the  beam  radius  r^  Is  smaller  than  the  conductivity  channel  r^,  the  beam  Is 

attracted  towards  the  channel  axis.  On  the  other  hand.  If  r^  >  r^  as  Is  usually 

true  In  the  beam  head,  we  have  shown  that  when  the  beam  resides  on  the  channel 

axis,  each  slice  will  see  an  average  force  repelling  It  from  the  channel.  But  In 

the  electrostatic  regime.  If  the  beam  Is  displaced  sufficiently  from  the  channel 

It  should  see  an  electric  dipole  attracting  It  back  toward  the  channel.  Using 

Dynasty  II  to  study  this  problem,  we  found  that  when  the  beam  displacement  Is  of 

the  order  of  the  beam  radius  and  r,  >  r  ,  a  stable  equilibrium  exists  with  the 

0  c 

beam  off  axis,  l.e.  Y  *  0.  The  equilibrium  displacement  from  the  channel  axis 
Increases  from  the  head  of  the  beam  to  its  tall.  A  detailed  description  of  these 
results  will  be  presented  In  a  later  report. 

When  an  external  electric  or  magnetic  force  F  acts  on  the  beam,  each  beam 
slice  will  be  subject  to  a  different  acceleration  F/my  if  y  is  a  function  of  i;. 


5 


If  the  shear  ia  F(?)/oy(5)  strong  compared  to  the  restoring  force  due  to  the 
beam's  self-pinch,  the  beam  will  tear.  The  break-up  will  continue  along  5  until 
for  some  value  the  restoring  force  becomes  stronger  than  the  sheared  external 
force;  from  this  point  on,  l.e.  ;  >  the  beam  will  hold  together.  Once  the 
"guiding  point"  is  determined,  we  can  find  (a)  the  deflection  due  to  the  force 
F,  (b)  what  portion  of  the  beam  is  torn  out,  and  (c)  the  value  of  i*e*, 

the  energy  of  that  part  of  the  beam  actually  guides  the  rest  of  the  beam. 

If  differs  substantially  from  7(5  *•  0),  then  a  problem  may  occur  in  aiming 

the  beam  subject  to  the  influence  of  this  external  force.  This  work  is  still  in 
progress  and  the  results  will  be  reported  in  a  later  report. 
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Appendix 

The  algorithm  for  inverting  a  bi-triagonal  system  with  scalar  coefficients 
is  well  knovm.^  We  shall  generalize  it  to  include  matrix  coefficients.  The 
equations  are 


for  1  <  i  <  R,  with  aj™^«  c_^®^“  0  for  1  <  m  <  4, 

1  R 


The  algorithm  is  as  follows: 


First  compute 


-  a^2),(3) 

®1  “i  *1  \-l  *1  \-l’ 


.O).  ,(3).  ,<3),a).  ,<.),(3,, 

e<‘>. 


with  0$”^“  b^“^  for  1  <  m  <  4, 
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